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Abstract 

Dauphas and Pourmand (2011) [Nature 473, 489-492] estimated the accre- 
tion timescale of Mars to be 1.8 l^'o Myr from the W isotopes of martian 
meteorites. This timescale was derived assuming perfect metal-silicate equi- 
libration between the impactor and the target's mantle. However, in the 
case of a small impactor most likely only a fraction of the target's mantle 
is involved in the equilibration, while only a small part of the impactor's 
core equilibrates in the case of a giant impact. We examined the effects of 
imperfect equilibration using results of high-resolution iV-body simulations 
for the oligarchic growth stage. These effects were found to be small as long 
as a planetary embryo has a deep liquid magma ocean during its accretion. 
The effect due to partial involvement of the target's mantle in equilibration 
is small due to the low metal-silicate partition coefficient for W suggested 
from the low Hf/W ratio of the martian mantle. The effect due to partial 
involvement of the impactor's core is also small because a large fraction of the 
embryo mass is delivered from small planetesimals, which are likely to fully 
equilibrate in the deep magma ocean on the embryo. The accretion timescale 
of Mars estimated by the Hf-W chronology is shorter than that expected for 
the minimum mass solar nebula model as long as more than 10% of each 
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impactor's core re-equilibrates with the martian mantle and the final stages 
of accretion are prolonged. This probably indicates that accretion of Mars 
proceeded rapidly due to solid and gas surface densities significantly larger 
than those for the minimum mass solar nebula or due to accretion of small 
fragments or pebbles. 

Keywords: Mars, impacts, Hf-W chronometer, iV-body simulations 



1. Introduction 

The theory of planet formation suggests that tens of Mars-sized embryos 
form from planetesimals in the inner solar system during the runaway and 
oligarchic growth stages, which last for ~ 0.1-10 Myr (Kokubo and Ida, 1998; 
Wetherill and Steward 1993). In the subsequent, much longer (~ 100 Myr) 
giant-impact stage, mutual collisions between embryos occur (Morishima et 
al., 2010; O'Brien et al., 2006). Mars is considered likely to be a remnant 
embryo because its mass is close to the theoretically predicted final mass of 
an oligarchic embryo, commonly referred to as the isolation mass (Lissauer 
et al., 1987), and its accretion timescale suggested from the Hf-W chronology 
(Dauphas and Pourmand, 2011; Nimmo and Kleine, 2007) is much shorter 
than that for the Earth-Moon system (Konig et al., 2011; Touboul et al., 
2007). The giant impact stage is rather stochastic, because the timing of 
the final large impact has a large dispersion in iV-body simulations starting 
with similar initial conditions (Morishima et al., 2010; O'Brien et al., 2006). 
However, the oligarchic growth stage is deterministic and the mass evolution 
of an embryo can even be expressed as a simple analytic solution if a uni- 
form size of planetesimals is adopted (Chambers, 2006). Therefore, if the 
accretion timescale of Mars is precisely determined by the Hf-W chronology, 
several important quantities which determine the accretion timescale may be 
retrieved, such as the gas and solid surface densities of the protosolar disk, 
and the planetesimal size. 

Dauphas and Pourmand (2011) estimated the accretion timescale of Mars, 
most precisely to date, to be 1.8 l^o Myr. This timescale was derived us- 
ing the isotopic evolution model of Jacobsen (2005) which assumes perfect 
metal-silicate equilibration between the impactor and the target's mantle 
during an impact (the more massive body is designated to be the target, 
while the less massive one is the impactor). However, two possible types 
of imperfect equilibration are conceivable; both of them result in a longer 
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accretion timescale. 

The first one is that only a part of the metallic core of the impactor is 
involved in equilibration. Mezger et al. (2012) showed that if ~ 40 % of each 
impactor's core re-equilibrates with the martian mantle, the estimated accre- 
tion time scale of Mars becomes a few times longer than in the case of perfect 
equilibration. Nimmo et al. (2010) investigated tungsten isotopic evolution 
of terrestrial planets using outputs of iV-body simulations of O'Brien et al. 
(2006) and found that in order to produce the terrestrial tungsten anomaly 
the degree of equilibration is required to be 30-80 %. Rudge et al. (2010) 
adopted various growth curves of the Earth and found that both perfect equi- 
librium and partial equilibrium model can reproduce the terrestrial tungsten 
anomaly and the abundances of siderophile elements. They also constrained 
the degree of equilibration to be at least 36%. 

Perfect equilibration requires emulsification of the impactor's core down 
to cm-scale droplets in the target's mantle (Rubie et al., 2003) because core 
formation is likely to occur even for a small impactor (down to ~ 10 km) 
due to heating by radiogenic decay of 26 Al and 60 Fe (Moskovitz and Gaidos, 
2011; Neumann et al., 2012). Whether sufficient emulsification occurs when 
a metallic core sinks in a liquid magma ocean was investigated using different 
types of hydrodynamical models (Dahl and Stevenson, 2010; Deguen et al., 
2011; Ichikawa et al., 2010; Samuel, 2012). Except for Dahl and Stevenson 
(2010), these studies consistently show that deformation and breakup of the 
iron core (neglected in Dahl and Stevenson's analytic model) is the main 
mechanism that leads to very efficient emulsification and metal-silicate equi- 
libration, as long as the core size is smaller than the depth of the magma 
ocean on the target. Kendall and Melosh (2012) showed that significant 
emulsification of the impactor's core already occurs during a high-velocity 
impact (they adopted an impact speed of 11.5 km/s). 

The second one, which has been paid less attention than the first one, 
is that only a part of the target's mantle is involved in equilibration. This 
is likely to occur if the impactor size is much smaller than the target size; 
such a situation is common during the oligarchic growth stage. Even if the 
impactor size is large, an upper metal-rich layer and a lower metal-poor layer 
may overturn due to the Rayleigh- Taylor instability, which may reduce the 
volume fraction of the target's mantle involved in equilibration (Sasaki and 
Abe, 2007). 

As discussed above, the impactor-to-target mass ratio is likely to be the 
important parameter for metal-silicate equilibration. The mass distribution 
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of impactors during the oligarchic growth stage is not clearly known. During 
the former runaway growth stage, the mass distribution can be described by 
a power-law with the exponent of ~ —2.5 (Kokubo and Ida, 2000; Ormel 
et al., 2010). During the oligarchic growth stage, large embryos separate 
from the continuous size distribution. In addition, mutual collisions between 
oligarchic bodies are expected. Chambers (2006) estimated that one-third of 
the mass of an oligarchic embryo is delivered by embryo-embryo collisions, 
provided that the mutual separation normalized by the Hill radius is fixed. 
This needs to be examined by direct iV-body simulations. 

In the present paper, we investigate accretion timescales of embryos and 
mass distribution of impactors, by conducting high resolution iV-body sim- 
ulations of the oligarchic growth stage. Some basics of the oligarchic growth 
stage are reviewed in Sec. 2. In Sec. 3, methods and results of two iV-body 
simulations are shown; one case with nebular gas and another case without 
gas. In Sec. 4, using the results of iV-body simulations, we model the isotopic 
evolution of the Hf-W system and examine the effects of possible imperfect 
metal-silicate equilibration. Discussion and summary are given in Sec. 5 and 
Sec. 6. 

2. Expected time evolution of the mass of an oligarchic body 

Before discussing the results of iV-body simulations, some basics of the 
oligarchic growth stage are briefly reviewed. Consider planetary embryos 
surrounded by small planetesimals. An embryo gravitationally influences 
planetesimals in an annulus around the embryo's orbit and these nearby 
planetesimals collide with the embryo. This annulus is called the feeding zone 
and its radial width is roughly equivalent to the orbital separation between 
neighboring embryos and is known to become 6 e r H as a result of orbital 
repulsion between embryos (Kokubo and Ida 1998), where b e ~ 10 is the 
scaling factor and r H is the Hill radius of the embryo given by 

ru = a ( ^V /3 = 2 i/3 ah (1) 
V 3M ©/ 

where a is the semimajor axis of the embryo, m c is the mass of the embryo, 
M & is the solar mass, and h is the reduced Hill radius used below. 

Defining the mass fraction of an embryo to the total solid mass in its 
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feeding zone to be / e , the mass of the embryo m e is given by 



_ 1/2 

m e = 2vr/ e a6 e r H S so i id = (2vr/ 6 S solid ) 3/2 a 3 ( 

-•«(t) S/2 (^) 3,2 (^) 3/2 (t^) 3/2 -» v 

Sgoiid is the solid surface density and mM ars is the mass of Mars (6.4185 xlO 26 
g). When an embryo sweeps up all solid material in its feeding zone (f e = 1), 
the embryo mass is called the isolation mass (Lissauer et al., 1987): 

m iso = m c (f c = 1). (3) 

The time evolution of f e is given by (Chambers, 2006) 

f = ^/ c 1/2 (l-/e), (4) 

where A is 

_ 31.7C / S solid \ 

where b c = 2 1 / 3 6 e , e is the orbital eccentricity of planetesimals normalized 
by h, P is the orbital period, p is the density of nebular gas. The factor C 
in Eq. §5§ represents acceleration of growth due to embryo-embryo collisions 
and we set C — 1.5 as estimated in Chambers (2006). If viscous stirring of 
embryos and gas drag are in equilibrium and the radius r of planetesimals is 
uniform, e is written as 



b c C D ap i 



where Cd is the drag coefficient assumed to be unity and p gas is the gas 
density. Eqs. ([5]) and mean that planetesimals frequently collide with an 
embryo, if e is low due to a large p gas or a small r or simply if S SO Hd is large. 
The solution of Eq. (j4j) is 



/ c (t) 1/2 = tanh 



'^dt' + atanh^O) 1 / 2 ) 



L^o 



(7) 
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and the time evolution of the embryo's mass is given as 



m e (t) = m- lso f c {tf/ 2 . (8) 
If A is independent of time and / e (0) is negligible, Eq. (JSJ) reduces to 



m e (t) = m; so tanh 3 



T. 



grow 



where 



T, 



grow 



2 

A' 



(9) 



(10) 



Eq. ([9]) shows that the embryo gains 44%, 90%, and 99% of its final mass 
during timescales lr grow , 2r grow , and 3r, 



' grow • 



The timescale r grow for Mars was 



estimated by Dauphas and Pourmand (2011) as 1.8 Myr by setting m iso 
to be the mass of Mars. The parameter A is usually time- dependent because 
gas dissipates with time. Even in such a case, an analytic expression can 
be obtained as long as A is analytically integrable with respect to time; for 
example, in the case where gas dissipates exponentially with time. 

To estimate / e (0), the embryo's mass at the transition from the runaway 
growth stage to the oligarchic growth stage is discussed. If the mass distri- 
bution is given by dn = km q dm (where n is the cumulative number of bodies 
inside the feeding zone of the embryo and k is a constant), q is ~ —2.5 at 
this transition (Kokubo and Ida, 2000; Ormel et al., 2010). The total mass 
inside the feeding zone is given by 



rriT 



q + 2 



m dn 

1 + 1 



mo 



-g-2 



m fi 



(11) 



where mo is the mass of the smallest planetesimal. In Eq. pip , we removed 
k using the definition of the largest body (embryo) given by 



00 , km q + l 
dn = — 

9 + 1 



(12) 



With m e ^> m and q = —2.5, the mass fraction of the embryo relative to 
the total mass is 



fu 



1 / m 



m T 3 V m. 



1/2 



(13) 
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Inserting Eq. (Ti"3"|) into Eq. 02]), the transition mass is given by 




6/7 



m tr = m e (f e = / tr ) = 



4/7 3/7 
iso m O 



(14) 



This mass coincides within a factor of two with the transition mass derived 
from the numerical simulations of Ormel et al. (2010) (their Eq. (13)). We 
adopt /c(0) = (m tr /m iso ) 2 / 3 from Eq. 

3. AT-body simulations 

3.1. Methods 

Two iV-body simulations were performed for systems in annuli around 
the current location of Mars, 1.5 AU; one with nebular gas (Sim. A) and the 
second one without gas (Sim. B). All input parameters are the same for both 
simulations except those for nebular gas. Each impact between two bodies 
is assumed to result in perfect merging. To keep the total surface density 
of solid bodies in a simulation annulus constant, the body supply boundary 
condition (Kokubo and Ida, 2000) is adopted; if the semimajor axis a of 
a planetesimal is larger than the outer boundary, a out , this planetesimal is 
sent to the inner boundary, a ini without changing its orbital eccentricity 
e and inclination i. With this boundary condition, planetesimals near the 
boundaries tend to confine embryos inside the annulus, mimicking viscous 
stirring of embryos outside the boundaries (Kokubo and Ida 2000). 

The total mass Md of planetesimals in the annulus is 2.0 x 10 27 g (~ 
3mMars)- The initial surface density of planetesimals is set to 



where E solid (1.5 AU) is adopted to be 5.0 g cm 2 . This is higher than that 
for the Minimum Mass Solar Nebula (MMSN; Hayashi, 1981) by 30%. The 
width of the simulation annulus Aa derived using the above Md and £ so iid is 
0.1896 AU. The boundary radii are then given as a- m = 1.5 AU — Aa/2 and 
a out = 1-5 AU +Aa/2. The initial number of planetesimals is 5,000 and the 
physical density p is 3.95 g cm -3 for all bodies. The initial mass distribution 
of planetesimals is given by a single power-law with q = —2.5. The initial 
mass ratio between the largest and the smallest bodies is set to be 20. This 
gives the mass of the smallest planetesimal m — 1.7 x 10 23 g. 




(15) 
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For the simulation with gas, the gas-to-solid ratio averaged over the annu- 
lus is set to be 240 at the beginning, as for the MMSN. The gas temperature 
is given as T = 220(r/1.5AU)- 1 / 2 K. The surface density of gaS, Zjg as , IS 
assumed to decay exponentially with a timescale of 2 Myr: 

/ a \ 7/4 / i \ 

s - (M) = (r^Au) exp (~^M^) Eg " (1 - 5AU ' 0) ' < 16 » 

where S gas (1.5 AU, 0) = 1,196 g cm -2 . The adopted decay time of the gas 
disk is consistent with observations (Fedele et al., 2010). In Eq. (fl6|) . the 
radial gradient of S gas is chosen so that the rotation velocity of the gas is 
exactly the local Kepler velocity at the midplane. Thus, there is no radial 
drift due to gas drag if e = i — 0. Aerodynamic gas drag and damping 
due to tidal interactions with the gas disk are taken into account following 
the approach of Morishima et al. (2010), whereas Type I migration (Gol- 
dreich and Tremaine, 1980) is neglected. A model of planetesimal formation 
in turbulence (Chambers 2010) suggests a planetesimal mass even smaller 
than m at 1.5 AU. To mimic smaller sizes of planetesimals in the gaseous 
disk, the drag force on the smallest planetesimal is enhanced by a factor 
of (mo/mci)^ 3 , where mci is the actual planetesimal mass of interest, and 
the enhancement factor smoothly decreases with increasing mass from itlq to 
10m (see Morishima et al., 2010). We set mci to be 10 21 g, which is close 
to that at 1.5 AU suggested by Chambers (2010). 

For this study, the iV-body code developed by Morishima et al. (2010) 
is used. In this code, the mutual gravity of all bodies is calculated with a 
parallel-tree method, and the orbital integration is calculated with a mixed- 
variable symplectic integrator. Since this code can handle a large number of 
particles and take a large time step, it is applicable to any stage of planetary 
accretion. The opening parameter for the tree method is 0.5 and the time 
step is 11.6 days. Each simulation was performed for a total duration of ~ 
14 Myr and took roughly two months using a single node with an eight-core 
processor on a supercomputer. 

3.2. Results 

Figures 1-3 show results for the simulation with nebular gas (Sim. A). 
Snapshots on the a-e plane are shown in Fig. 1. The large embryos are 
displayed as red filled circles with horizontal branches with a length of 5rn 
on each side. We define an embryo as a body more massive than 10 25 g 
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(~ m tr (Eq. (j!4p ) for b e = 10), but bodies less massive than one-tenth of the 
most massive body's mass are excluded, as such a small body behaves like 
a planetesimal under gravitational stirring by massive embryos. As embryos 
grow with time, the number of planetesimals decreases and their eccentricities 
increase. At t — 1 Myr, many orbital overlaps of embryos can be observed. 
Through collisions between embryos, two large embryos form. Their orbital 
separation is roughly IO/h as found in Kokubo and Ida (1998, 2000). At the 
end, their masses reach 80-85% of the mass of Mars. 

The cumulative number of bodies and the velocity distribution as a func- 
tion of mass are shown in Fig. 2. While the power-law index for the mass 
distribution is initially —2.5, it gradually evolves to ~ —2 as a result of em- 
bryo growth. The large embryos have low e due to tidal interactions with 
gas and dynamical friction of surrounding planetesimals. For small planetes- 
imals, e decreases with decreasing m due to gas drag. The analytic estimates 
of the velocity dispersion (Eq. (jSJ)) of the smallest bodies are shown on the left 
side as dotted lines and coincide well with those of the A-body simulation. 

The masses of the two surviving embryos as a function of time are shown 
in the top panel of Fig. 3. In the same panel, the impactor-to-embryo mass 
ratios (< 1) are shown as diamonds for embryo-embryo collisions. The total 
mass gained by embryo-embryo collisions, / e -e, relative to the mass of the 
embryo is shown in the middle panel of Fig. 3 for each embryo. The number 
of embryos N e and the mean separation normalized by the Hill radius b e are 
shown in the bottom panel of Fig. 3. Initially, N e increases and b e decreases 
with time. After b e reduces to ~ 5 around 1 Myr, embryo-embryo collisions 
start to occur frequently and N e decreases while b e increases. As a result, 
surviving embryos gain roughly half of their masses by embryo-embryo col- 
lisions. These fractions decrease with time near the end of the simulation 
because all other embryos are swept up at 6 Myr and thereafter only small 
planetesimals collide with the embryos. The fractions extrapolated to the 
end state where all planetesimals are completely swept up are 0.1-0.2. 

On the top panel of Fig. 3, we also plot the theoretically expected mass 
evolution (Eq. OH])). The velocity dispersion calculated from Eq. §6§ increases 
with time and its upper limit is fixed to be the escape velocity of embryos 
in Eq. ©. We adopt b c = 9.1 in Eq. © so that m iso = M d /2 = 10 27 g; this 
is the expected final mass of the two surviving embryos. These conditions 
give r grow = 3.06 Myr in the beginning of the simulation. The iV-body 
simulation shows a good agreement with the analytic estimate up to ~ 3 
Myr. Good agreement with the analytic model adopting a bimodal size 
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distribution is reasonable as a large mass fraction is in smallest planetesimals 
in the simulation. 

However, after ~ 3 Myr the simulation shows that the embryos grow by a 
factor of three slower than the analytic estimate. There are several plausible 
reasons reducing the growth rate in the simulation. First, embryo-embryo 
collisions cease around 3 Myr while they are assumed to occur constantly in 
the analytical estimation (C is fixed to be 1.5 in Eq. ([5])). Cease of embryo- 
embryo collisions in the middle of planetary accretion is also seen in N- 
body simulations without the periodic boundary condition (Morishima et 
al., 2008). Second, the surface densities in the feeding zones of embryos are 
reduced relative to those near the boundaries of the annulus, due to strong 
scattering by embryos. A non-uniform spatial distribution of planetesimals is 
seen even without gas (Ida and Makino, 1993) but non-uniformity is enhanced 
with gas drag (Tanaka et al., 1997). This effect is already seen at t — 2 Myr 
in Fig. 1 and also seen in simulations of Kokubo and Ida (2000). Third, the 
width of the feeding zone of an embryo increases with orbital eccentricities 
of planetesimals. Thus, during the late stage, a non-negligible part of the 
feeding zone lies outside the simulated annulus. This effect is likely to be 
relatively unimportant compared to the first and second effects. Overall, the 
analytic model is likely to underestimate the timescale of sweeping up of 
remnant planetesimals during the late oligarchic growth stage. 

Figures 4 and 5 show the same as Figs. 1 and 3 but for the simulation 
without nebular gas (Sim. B). Overall evolution of Sim. B is similar to Sim. A. 
Due to a higher velocity dispersion in the absence of gas, however, embryo 
growth is significantly slower than in Sim. A. The velocity dispersion is found 
to be roughly half of the escape velocity of the largest embryo throughout the 
simulation. The final number of embryos is only one in Sim. B in contrast 
to two in Sim. A. This is because orbits of embryos become more unstable 
without damping due to tidal interactions with the gas disk (Iwasaki et al., 
2002; Kominami et al., 2002). It is expected that the embryo-embryo spacing 
becomes much larger without gas than with gas in the late stage if the simula- 
tion annulus is much wider. Related to this, the peak number of embryos for 
Sim. B (Fig. 5) is only half of that for Sim. A. Despite these differences, the 
mass fraction of the largest embryo obtained by embryo-embryo collisions in 
Sim. B is similar to those in Sim. A. The mass growth curve roughly coincides 
with the analytic estimation assuming that the velocity dispersion is half of 
the escape velocity of the largest embryo (Fig. 5) and m iso = M d = 2 x 10 27 
g (b e = 14.4 and r grow = 15.63 Myr). The growth rate in the late stage of 
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the simulation is likely to be lower than the analytic estimate for the same 
reasons mentioned for Sim. A. 

4. Evolution of Hf-W isotopes 

The Hf-W chronology provides strong constraints on core formation tim- 
ing in planetary bodies (Jacobsen, 2005; Kleine et al., 2009). 182 Hf decays 
to 182 W with a half-life of 8.9 ±0.1 Myr. Hf and W are both refractory 
elements and their relative abundances in bulk planetary bodies should be 
chondritic. During core formation, lithophile Hf entirely remains in a silicate 
mantle whereas siderophile W is preferentially partitioned into a metallic 
core. 

The evolution of Hf-W isotopes is investigated using the results of the N- 
body simulations. Since the iV-body simulations were halted before complete 
sweep up of planetesimals, we extrapolate the growth curves of the embryos 
(dashed curves in Figs. 3 and 5) using Eq. (J4J); the growth timescale 2/ A of 
an embryo is derived using the masses of the embryo at t — 10 Myr and at 
the end of the simulation (t ~ 14 Myr) assuming A is constant during this 
period. As an extreme comparison, we also calculate the evolution of the 
Hf-W isotopes for without subsequent accretion after the end of the 

simulation assuming that remaining planetesimals are suddenly dispersed. 

4-1. Formulation 

We make the following assumptions: 

1. The abundances of Hf and W in all bulk bodies are chondritic. 

2. The metal-silicate partition coefficient for W, _D W , is independent of 
time and space. The silicate fraction y and the complementary metal 
fraction 1 — y are also constant for all bodies. 

3. All bodies experience core formation at the beginning of the solar sys- 
tem due to 26 Al radiogenic heating and their cores and mantles are fully 
equilibrated at this time (the timing of the core-mantle differentiation 
is found to have little effect). 

4. During a collision, the core of the target (the more massive body) is 
not involved in metal-silicate equilibration. The target's mantle and the 
impactor's core are partially involved in equilibration and the fractions 
of involved masses are defined to be k tm and k- lC) respectively. 
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MO = ( ,^"L," ,W u, - 1 ) x 10*. (17) 



5. After equilibration, the impactor's core merges with the target's core. 
The target's mantle is well mixed and homogenized immediately after 
the core merging, as well as the core, but there is no material exchange 
between the core and the mantle. 

The W isotope ratio 182 W/ 183 W for a sample relative to that for the 
CHondritic Uniform Reservoir (CHUR) is defined as 

> ^W/ 183 W)(t) 

W/ 18 W) C hur(*) 

where t is the time measured from the beginning of the solar system (4,568 
Myr ago; Kleine et al., 2009). The value for the present-day Mars is estimated 
from Shergottites to be Ae = 2.68 ± 0.26 (Dauphas and Pourmand, 2011). 

The radiogenic change of Ae during a time interval between two colli- 
sions (the first and the second collisions occurring at tl and £2) is given by 
(Jacobsen 2005) 

Ae(t2) - Ae(f 1) = C w {e~ xtl - e- xt2 ) , (18) 
where A is the decay constant of 182 Hf. The coefficient C w is given as 



182TTf \ *=° 

nL J CHUR 



where 

/18( 

n „ T = 1 n 4 I — 

2 W 



180 j^f \ t=0 



= 10 ( ISw ) ( 20 ) 



CHUR 



f „,/w = ( rm/^w) _ A (21) 

1 V( 18 "Hf/ 183 W)cHUR )' 1 ' 

Strictly speaking, g w should be the value at time t but can be well ap- 
proximated as a constant for the Hf-W system. The values from Dauphas 
and Pourmand (2011; and references therein) are used: gw — 1-07 x 10 4 , 

(182 H f/18O Hf )t=0 uR = q 72 x 10 -5 > jHf/W = 3 33 ± Q 5Q ^ &nd A = Q.0779 

Myr -1 . Only the uncertainty in f m / w is taken into account, as it has the 
largest effect on the estimated accretion timescale. With assumptions 2 and 
3, f m / w is time-independent for all bodies. Additionally using assumption 
1, / Hf / W is given as (Jacobsen, 2005) 

yHf/w _ (1 — y)D w 
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A change in Ae due to a collision, ignoring radiogenic decay during the 
collision, is 

Aei = / t Ae + fiAe u (23) 



where Aeo and Aei are Ae of the target before and after the collision and Ae; 
is Ae of the impactor. The reduction factor f t represents how much 182 W 
in the target's mantle is transported to the target's core by the impactor 's 
core while /; represents how much the super chondritic impactor's mantle is 
added to the target's mantle. These factors are given as (see Appendix) 

- = (fct m + ff) + (i-MMJ Hf/w ff , , 

h (k tTa + g + k ic f^g)(l + g) 1 ) 

and 

, = (1 - k ic )(k tm + g)g 

Ji (h m + g + k ic f^g)(l + g y 1 ) 

where g(< 1) is the mass ratio of the impactor to the target. In the case 
of k ic = 1, fi will become zero. Thus, a change in Ae does not depend 
on the isotopic fractionation history of the impactor. For k lc = k tm = 1, 
Eq. fH]) is reduced to Eq. (A6) of Kleine et al. (2009). This case is called 
the mantle equilibration scenario in Nimmo and Agnor (2006). The core- 
merging scenario of Nimmo and Agnor (2006) corresponds to with 
k [c = and k tm = 1. 

Time evolution of Ae of an embryo is calculated using Eqs. (fl8l) and ( 1231) . 
After the end of an iV-body simulation, impacts to the embryo are assumed 
to occur every At and the impactor mass is given by m e (t) —m e (t — At). We 
adopt a somewhat large At(= 10 6 yr) to represent remaining medium-sized 
impactors (g ~ 0.01). For k ic < 1, Ae's of all bodies are calculated while we 
assume that Aei = Aeo after the end of the simulation. 

4-2. Results 

We first show the isotopic evolution for the case of perfect equilibration 
{ktm = k ic = 1). Next, the effect of partial equilibration of the target's 
mantle is examined (k tm < 1). Third, the effect of partial equilibration of 
the impactor's core is studied (k[ C < 1). Finally, partial equilibration of both 
the target's mantle and the impactor's core is taken into account. 
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4-2.1. Perfect equilibration 

Figure 6 shows the time evolution of Ae for the case of perfect equilibra- 
tion (ktm = k- lc = 1). Results for two embryos from simulation A (red and 
blue solid curves) and for one embryo from simulation B (black curve) are 
shown. The dashed curves show a case without accretion after the end of the 
iV-body simulation. With increasing time, Ae increases due to radiogenic 
decay of 182 Hf to 182 W, while impacts reduce Ae as 182 W is transported from 
the mantle to the core. The larger the impactor-to-target ratio g, the larger 
the decrease in Ae. After ~ 20 Myr, Ae starts to decrease because decrease 
due to impacts exceeds increase due to radiogenic decay. The present-day 
values of Ae become less than 0.3 for all cases and are much lower than the 
martian value. If accretion is halted at the end of the iV-body simulations, 
the present-day values of Ae become much larger but still lower than the 
martian value. 

4-2.2. Partial equilibration of the target's mantle 

If the impactor size is much smaller than the target size, the entire target's 
mantle may not be involved in the metal-silicate equilibration. This effect 
is modeled by assuming that the volume of a portion of the target's mantle 
interacting with the impactor's core is proportional to the impactor's volume: 

k tm = mm[c tm g, 1.0]. (26) 

where c tm is the proportionality coefficient (ignoring a factor of y). A similar 
form is also used in Sasaki and Abe (2007). If the target has a liquid magma 
ocean on the surface, the volume of the interacting portion is probably given 
by the product of the crater area and the magma ocean depth. In this case, 
c tm is likely to be much larger than unity particularly for small impactors. 
On the other hand, if the target surface temperature is lower than the solidus 
temperature of the rocks composing the mantle, the impact-induced melt in 
the target's mantle may be taken as the interacting portion. The volume 
of the impact-induced melt primarily depends on impact velocity, and if the 
impact velocity is about the escape velocity of Mars, only the isobaric core, 
whose volume is a few times of the impactor volume (c tm ~ 3), results in 
melting (Tonks and Melosh, 1993, 1992). 

Figure 7 shows the present-day value of Ae as a function of Ct m . The effect 
of partial equilibration is found to be unimportant for large c tm . Inserting 
fc tm = c tm g into Eq. (1211) and adopting c tm ^> / Hf / W ; we have f t ~ (1 — 
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f m ' w g)/(l + g). This is equivalent to the case of perfect equilibration if 
g <C 1. This indicates that the tungsten concentration in the impactor's core 
tends to saturate if the volume of the equilibrating silicate portion exceeds 
the criterion, Ct m ~ J Hf / w . The present-day value of Ae becomes much larger 
if c tm < / Hf / W (ft ~ (1 - c tm g)/(l + g) for c tm « / Hf / W ). Such a case is 
unlikely for Mars because its / Hf / W is low and c tm is likely to be at least ~ 3 
as discussed above. The effect of partial equilibration may be much more 
important for the Earth which has a much larger J Hf / w ; the Hf/W ratio of 
25.8 (Kdnig et al., 2011) gives / Hf / W = 25.3 using the Hf/W ratio of 0.98 for 
CI chondrites. 

4-2.3. Partial equilibration of the impactor's core 

Cases of a constant k ic (< 1) for all impacts are first considered while 
keeping k tm = 1. The reductions of Ae after impacts are much smaller 
than in the case of perfect equilibration. Eq. (T24]) shows that f t = (1 + 
9 + k\cf m ^ g)~ l ■ Thus, the effect of partial equilibration of the impactor's 
core is equivalent to reducing f m / w to k ic f m ^ w , or to reducing the metal- 
silicate partition coefficient D w to k ic D w . In addition, the disequilibrated 
impactor's mantle is added to the target's mantle (Eq. ( 12"5"|) ). This leads to 
an additional increase in Ae. 

Figure 8 shows the present-day value of Ae as a function of k[ C . As 
discussed, Ae increases with decreasing fc ic . In the limit of k lc — > 0, Ae = Cw 
(Eq. ( lT9"j) ) regardless of accretion history. It is found that Ae for present-day 
Mars is fulfilled with k- lc less than 0.1. This value is even smaller than the 
lower limit (0.36) suggested for the Earth (Rudge et al., 2010). If accretion 
is halted at the end of the iV-body simulation with gas, k lc needs to range 
from 0.2 to 0.7. 

4-2.4- A model with a deep magma ocean 

In this section, we consider where Mars has a deep magma ocean 

during its accretion. Whether this assumption is reasonable is discussed in 
the next section. As found in hydrodynamic modeling (Deguen et al., 2011; 
Samuel, 2012), ki c is likely to decrease with increasing impactor size. As a 
nominal model, we adopt the following form: 

h c = minfcicgf" 1 , 1], (27) 

where c lc is the proportional coefficient and we adopt c lc = 0.01. This means 
that large impactors only partially equilibrate with the target whereas smaller 
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impactors (g < 0.01) fully equilibrate. This model is probably close to the 
condition derived by Samuel (2012), who found that impactors with sizes 
smaller than the depth of the target's magma ocean fully equilibrate. We 
also adopt partial equilibration of the target's mantle using Eq. fT26|) with 
Ctm = 10. As discussed in Sec. 4.2.2, this is likely to cause only a very small 
difference from the case of /c tm = 1. 

The evolution curves of Ae are shown in Fig. 9. There is no large decrease 
in Ae at impacts unlike in Fig. 6 because k- lc for large impactors is very 
small. Nevertheless, the increase in Ae as compared with those for the cases 
of perfect equilibration is rather small, because more than half of an embryo 
mass is delivered by small impactors and their cores fully equilibrate with 
the embryo mantle. The present-day values of Ae are close to the values 
obtained with k\ c = 0.5 — 0.7 in Fig. 8. This is reasonable as the fraction 
of small impactors is roughly two-thirds as shown in Figs. 3 and 5. In any 
cases, the present-day values of Ae are much lower than the martian value. 
If accretion is halted at the end of the iV-body simulation with gas, the 
present-day values of Ae becomes barely consistent with the martian value. 

So far, we have fixed the value of f m / w to be 3.38. Nimmo and Kleine 
(2007) showed that different values of f m / w result in very different accretion 
timescales. Thus, the uncertainty of f m / w needs to be taken into account, 
although it is now much smaller owing to the work of Dauphas and Pourmand 
(2011). Figure 10 shows Ae as a function of J Hf / w . It is found that Ae is 
almost independent of f m / w (solid lines). Increase of Ae with time due to 
radiogenic decay is simply proportional to f m / w (Eq. (TIB"]) ). On the other 
hand, the decrease of Ae due to an impact also increases with increasing 
yHf/w ^ see f[24l) and Sec. 4.2.2). If the accretion timescale is longer than 
the radiogenic decay timescale, both effects roughly compensate each other 
and the present-day Ae depends very weakly on f m / w . On the other hand, if 
accretion rapidly completes as in the cases shown by dashed lines in Fig. 10, 
the present-day Ae is primarily determined by radiogenic decay so Ae is 
roughly proportional to f m / w . Similar dependence and independence can 
be found if f m / w is varied in calculations shown in Fig. 6. 

5. Discussion 

In Sec. 4.2.4, we showed that as long as a Mars-analog has a deep magma 
ocean during its accretion and the final stage of accretion is prolonged, the 
present-day Ae turns out to be much lower than the martian value. This 
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may indicate that accretion of Mars proceeded and terminated much more 
rapidly than in our model calculations, because (1) the solid surface density 
was larger, (2) the gas surface density was larger, or (3) the planetesimals 
size was smaller than what we have assumed. For the first possibility, the 
terrestrial planets might have accreted in a narrow annulus around 1 AU, 
and Mars is ejected to the current location (Hansen, 2009; Morishima et al., 
2008; Walsh et al., 2011). Although this scenario is favorable for the strong 
radial mass concentration of the terrestrial planets, the Earth is likely to 
accrete too rapidly to explain its Ae and it may be difficult to produce the 
large difference in the oxidization states (represented by / Hf / W and the FeO 
content) between Earth and Mars (but see Morbidelli and Rubie (2012)). 
For the second possibility, the high gas density may be achieved in a low 
turbulent viscosity region called the dead zone (Gammie, 1996; Morishima, 
2012). For the third possibility, growth of embryos is significantly accelerated 
due to accretion of fragments of planetesimals (Chambers, 2006; Kobayashi 
et al., 2010; Wetherill and Stewart, 1993) or small pebbles that are precursors 
of planetesimals (Lambrechts and Johansen, 2012; Morbidelli and Nesvorny, 
2012). A quantitative discussion assuming perfect metal-silicate equilibration 
can be found in Kobayashi and Dauphas (2012). The mechanisms listed 
above also help to terminate accretion quickly because Mars is ejected to a 
low surface density region, or remnant planetesimals are removed by strong 
gas drag. Perturbations from giant planets also help to disperse remnant 
planetesimals (e.g., Morishima et al., 2010). Unfortunately, these effects 
cannot be properly evaluated in local simulations with the periodic boundary 
condition, although it is evident that rapid dispersal of remnant planetesimals 
significantly enhances tungsten anomaly (dashed lines in Fig. 6 and 9). 

Alternative to the rapid accretion discussed above, the degree of metal- 
silicate equilibration might have been very low (k ic < 0.1, see Sec. 4.2.3) 
while accretion of Mars was as slow as in our simulations (see also Mezger et 
al., 2012). Such a low degree of equilibration might be feasible if accretion 
still proceeds while the global magma ocean starts to solidify. However, 
geochemical evidence seems to suggest that Mars gained most of its mass 
while it had a deep magma ocean. The most recent study of partitioning 
of the siderophile elements indicates a deep (> 1,000 km) global magma 
ocean on Mars (Righter and Chabot, 2011), contrary to a shallow magma 
ocean suggested by previous studies (Righter and Drake, 1996; Kong et al., 
1999). A deep magma ocean of Mars is also indicated from a differentiation 
model (Elkins-Tanton et al., 2005) that can produce magma source regions 
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consistent with element data for SNC meteorites and the martian crust. For 
Mars-size bodies, both impact heating and radiogenic heating due to the 
decay of 26 Al are likely to contribute to formation of magma oceans (Rubie 
et al., 2007; Sramek et al., 2012). Numerical models suggest that the presence 
of a magma ocean on Mars is very likely (Monteux et al., 2010; Ricard et al., 
2009; Sramek et al., 2012) but its depth and duration depends on various 
uncertain parameters. An important parameter is impactor sizes which can 
be directly obtained from iV-body simulations. Thermal evolution of the 
embryos in our iV-body simulations is investigated in a companion paper 
(Golabek et al., in preparation). 

A late stage giant impact is favorable for the origin of the geological 
dichotomy of Mars (Golabek et al., 2011; Marinova et al., 2008; Nimmo 
et al., 2008; Wilhelms and Squyres, 1984). The impact is unlikely to have 
significantly reset the chronometer probably because the impactor size was 
not large enough to stir the entire martian mantle or emulsification of the 
impactor 's core was inefficient. 

6. Summary 

We have performed high-resolution iV-body simulations for oligarchic 
growth of Mars, using the body supply boundary condition which keeps the 
total solid mass in the simulation annulus constant. One simulation in the 
gas-free environment and another simulation with nebular gas were carried 
out. The surface densities of solid and gas are larger than those for the 
minimum mass solar nebula model by only 30%. The time evolution of the 
embryo mass in nebular gas coincides well with that of the analytic model of 
Chambers (2006) until the middle of accretion. However, the iV-body simu- 
lations show much slower accretion than the analytic model in the late stage. 
The mass delivered by embryo-embryo collisions relative to the embryo mass 
is about half at maximum in the middle of accretion. This fraction decreases 
to 0.1-0.2 when planetesimals are completely swept up because only small 
planetesimals collide with the embryos during the late stage of accretion. 

We have calculated the Hf-W isotopic evolution of the embryos, using 
accretion histories from the iV-body simulations. It is likely that only a part 
of the target's mantle is involved in the equilibration if the impactor size is 
small while only a small fraction of the impactor's core is involved for a large 
impactor. These effects are modeled and examined using the output of the 
A-body simulations. It was found that as long as the target has a deep liquid 
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magma ocean during its accretion, the effects of imperfect equilibration are 
small. The effect of partial involvement of the target's mantle is small as 
long as the volume of the target's mantle portion involved in equilibration 
relative to the impactor volume is larger than f m / w . This is likely to be the 
case for Mars as its f m / w value is low. The effect of partial involvement of 
the impactor's core is also small because a large fraction of the embryo mass 
is delivered by small impactors, which are likely to fully equilibrate in the 
deep magma ocean of the embryo (Samuel, 2012). 

We found that with a prolonged accretion suggested from our simulations 
(~ 100 Myr), the high value of the martian tungsten anomaly can be achieved 
only if less than 10% of each impactor's core re-equilibrates with the martian 
mantle. This indicates that growth of Mars proceeded rapidly due to solid 
and gas surface densities significantly larger than those for the minimum 
mass solar nebula or due to accretion of small fragments or pebbles. 
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Appendix: changes of tungsten isotope concentrations during a 
collision 

The mass concentrations of the radiogenic and stable isotopes of tungsten 
are defined as 182 W and 183 W. Let us assume that a change in 182 W (or 183 W) 
due to a collision is given as: 

182 W tmil = A 182 W tm , + /i 182 W im + / c 182 Wchur, (28) 

where / t , /i, and f c are the coefficients derived below, the subscripts tm, 
im, and CHUR represent the target's mantle, the impactor's mantle, and 
the CHondritic Uniform Reservoir, and the subscript numbers and 1 for 
the target's mantle mean before and after the collision, respectively. With 
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assumptions 2 and 3 from Sec. 4.1, the concentration of the stable isotope 
183 W does not change due to a collision even for k[ C < 1 or k tm < 1 (proved 
below). This leads to 

183 W W = 183 W tm , = 183 W im = - f f c / 83 W C hur, (29) 

1 — Jt — Ji 

and 

Ae a = / t Ae Q + fiAe h (30) 

where Aei and Aeo are Ae of the target's mantle after and before the collision 
and Ae; is for the impactor's mantle. 

In the following, the factors f t , and f c are derived. The total mass 
involved in equilibration is 

M eq = k tm y t M t + y-M + k ic (l - Vi )M u (31) 

where M t and M ; are the masses of the target and the impactor, and y t and 
yi are the silicate fractions of the target and the impactor. Since the bulk 
composition of the impactor is assumed to be chondritic, we have 

182 W C hur = 2/i 182 W im + (1 - yi) 182 W ic , (32) 

where the subscript ic denotes the impactor's core. Using Eq. (|32l . the total 
mass of equilibrated 182 W supplied by the impactor relative to the impactor 
mass is given as 

182 w iiCq = yi 182 w im + Mi-yi) 182 w ic 

= (1 - k ic ) yi 182 W im + fc ic 182 W C HUR. (33) 

The concentration 182 W eqi o averaged over the masses involved in the equili- 
bration is given by 

182 W eq , = ^(fc tm2 / t M t 182 W tm ,o + Mi 182 W i)eq ) (34) 

The silicate fraction of the equilibrated masses is 

k tm y t M t + yiMj 

*> = K, (35) 
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After the equilibration, the concentration 182 W eqi i of the equilibrated silicate 
portion is (Nimmo and Agnor, 2006) 

182 W 

182 W _ VV cq.O (oc\ 

The fraction of 182 W tmi i for the entire target's mantle is given by averaging 
the equilibrated and non-equilibrated portions as 

i82 W _ hmytM t + y-M 182 (1 - ktm)ytM t m ( 

W tm,l - T-j—. 7^— W eq ,l H r-j— — W t m,0 ' 

y t M t + yM\ VtM t + yMx 

From Eqs. f[3"4"]) -f l3"T]) . we have 

= iM)c, M -" (40) 

where 

Cj = k tm y t M t + y x M h (41) 

C 2 = k ic D w (l- yi )M h (42) 

C 3 = y t M t + y-M-, (43) 

For y = & = y t , Eqs. §E§ and ([39]) are reduced to Eqs. ([21]) and fl25]) . 
respectively, using Eq. (j22J). 

Eq. (j3"7j) is also applicable to the stable isotope 183 W. Consider the case 
where both the impactor and the target body have not experienced any 
collision but have experienced core formation with full metal-silicate equili- 
bration. For these bodies with y = y\ = y t , we have 

183 w 

183 W _ 183 W W C HUR , ... 

W tm , - W im - y + (i _ y)DW . (44) 

Inserting this equation into Eq. ([37]) . we obtain 183 W tm ,i = 183 W tm , (Eq. ([29]) ). 
This relationship holds for subsequent collisions. Therefore, the concentra- 
tion of the stable isotope is not changed by collisions, as long as D w and y 
are constant. 
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Figure 1: Snapshots in the a-e plane for the iV-body simulation with nebular gas (simula- 
tion A). The circles are proportional to the radii of planetesimals and embryos. Embryos 
are displayed as red filled circles with horizontal branches with a side length corresponding 
to 5rn. The total numbers of bodies are 3,397, 2,721, 2,012, and 954 from the top panel 
to the bottom. 
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Figure 2: Cumulative number and orbital eccentricity as a function of mass for simula- 
tion A. The distributions at the four times from Fig. 1 are plotted. The dashed and dotted 
lines in the upper panel are slopes with q = —2.0 and —2.5. The dashed lines in the lower 
panel show the analytically estimated values for the smallest planetesimal. 
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Figure 3: (Top) The time evolution of the masses of the two surviving embryos relative to 
the isolation mass (10 27 g) for simulation A. Diamonds indicate the impactor-to-embryo 
mass ratios (< 1) for embryo-embryo collisions. The dotted curve is the embryo mass 
derived from the analytic model. The vertical dotted line indicates the time at the end 
of simulation. In the top and middle panels, the red (blue) color is used for the inner 
(outer) surviving embryo, and the extrapolated curves until the end of compete sweep 
up of planetesimals are shown by dashed curves. (Middle) Mass fractions of the embryos 
obtained by embryo-embryo collisions. (Bottom) Number of embryos N e (black line) and 
the mean separation b c (blue line) of neighboring embryos normalized by the mutual Hill 
radius. 
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Figure 4: Same as Fig. 1, but for the simulation without gas (simulation B). There is only 
one surviving embryo in this case. The total numbers of bodies are 3,415, 2,874, 2,333, 
and 1,392 from the top panel to the bottom. 
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Figure 5: Same as Fig. 3, but for simulation B. The isolation mass is set to be 2 x 10 27 
g. For the calculation of the analytic growth curve (dotted curve) in the top panel, the 
velocity dispersion is assumed to be half the escape velocity of the largest embryo. 
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Figure 6: Time evolution of Ae for the surviving embryos for the case of perfect equilibra- 
tion between the target's mantle and impactors (fc tm = k- lc = 1). The red and blue curves 
are for the inner and outer embryos from simulation A whereas the black curve is for 
the embryo from simulation B. The solid curves are cases where the growth curves of the 
embryos are extrapolated to complete sweep up of planetesimals while the dashed curves 
are cases where accretion is stopped at the end of TV-body simulations (t ~ 14 Myr) . The 
error bar at the upper right represents Ae for present-day Mars. 
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Figure 7: Present-day values of Ae for the surviving embryos for the case with partial 
equilibration of the target's mantle (/c tm < 1.0 and k lc = 1.0). Eq. (|26|) is used for 
modeling of k tul . The meanings of curves are the same as in Fig. 6. The region between 
the two horizontal dashed lines represents the range of Ae for present-day Mars while the 
vertical dashed line represents ct m = / Hf / W = 3.38. 
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Figure 9: Same as Fig. 6 but for the case with partial equilibration of both the target's 
mantle and the impactors' cores. Eq. (|26l) with c tm = 10 is used for the modeling of fc tm 
whereas Eq. (|27|) with c; c = 0.01 is used for the modeling of k lc . 
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Figure 10: Present-day values of Ae for the surviving embryos as a function of / Hf / W 
The partial equilibration model from Fig. 9 is used. The meanings of the solid and dashed 
lines, and the horizontal dotted lines are the same as those in Fig. 7. The region between 
the two vertical dotted lines represents the range of f m / w estimated by Dauphas and 
Pourmand (2011). 
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